#!/bin/sh
#
# This is the script for generating files for a specific Dalton test job.
#
# For the .check file ksh or bash is preferred, otherwise use sh
# (and hope it is not the old Bourne shell, which will not work)
#
if [ -x /bin/ksh ]; then
   CHECK_SHELL='#!/bin/ksh'
elif [ -x /bin/bash ]; then
   CHECK_SHELL='#!/bin/bash'
else
   CHECK_SHELL='#!/bin/sh'
fi

#######################################################################
#  TEST DESCRIPTION
#######################################################################
cat > ccsdmm_modelmix.info <<%EOF%
   ccsdmm_modelmix
   ---------------
   Molecule:         H2O
   Wave Function:    CCSD / 3-21G
   Test Purpose:     Check energy, dipole moment quadrupole moment,
                     first ground state excitation energy, transition
                     moment and static polarizability calculated
                     using the SPC/SPC_EC1/SPC_EC3/SPC_E01 models for
                     the 12 water molecules surrounding the QM
                     water molecule. OLDTG = .TRUE.
%EOF%

#######################################################################
#  MOLECULE INPUT
#######################################################################
cat > ccsdmm_modelmix.mol <<%EOF%
ATOMBASIS
QM/MM H2O(QM)+ 5 H2O(MM)
------------------------
    4    0         1 1.00D-15
        8.0   1    Bas=3-21G
O            0.000000        0.000000        0.000000 0 1
        1.0   2    Bas=3-21G
H           -0.756799        0.000000        0.586007 0 2
H            0.756799        0.000000        0.586007 0 3
   -0.669    12    Bas=MM
O           -6.022295       -6.249876       -2.389355 1 1  # CORD(I,J), ISUBSY(I), ISUBSI(I)
O           -0.590747        4.825666       -1.709744 2 1  # coordinate, sub system, sub system site
O            2.365069       -0.266593        1.169946 3 1
O           -8.979615       -1.935917       -5.707554 4 1
O           -5.696915       -2.203270        0.274131 5 1  #Subsystem 5, subsystem site 1
O           -3.468351       -4.870005       -5.943324 6 1
O           -7.496926       -3.470414       -3.508220 7 1
O           -0.535706       -3.459430      -10.918135 8 1
O           -4.219513        2.293357        4.687821 9 1
O            0.893082       -5.976085       -1.522108 10 1
O            0.577780       -0.898605       -9.957619 11 1
O           -7.330474       -0.800736       -1.814272 12 1
   0.3345    24    Bas=MM
H           -6.934736       -6.264225       -2.100595 1 2
H           -5.562128       -6.801680       -1.756974 1 3
H           -1.493716        5.142033       -1.683118 2 2
H           -0.065506        5.576575       -1.433333 2 3
H            2.261479       -1.103382        1.622925 3 2
H            3.132239       -0.390013        0.611059 3 3
H           -8.312439       -2.243695       -5.094117 4 2
H           -8.719849       -2.307644       -6.550450 4 3
H           -4.847884       -2.109424       -0.157695 5 2 #Subsystem 5, subsystem site 2
H           -6.336813       -2.008951       -0.410639 5 3 #Subsystem 5, subsystem site 3
H           -4.314395       -5.108436       -5.564490 6 2 #-----------------------------
H           -3.655119       -4.101117       -6.481904 6 3 # Together with subsystem site 1 
H           -7.274548       -4.398536       -3.580814 7 2 # they form a water molecule,
H           -8.166766       -3.437216       -2.825307 7 3 # which is labeled subsystem 5!!!
H            0.271775       -3.713098      -10.471183 8 2
H           -1.195126       -3.432889      -10.224879 8 3
H           -3.911733        1.437798        4.986865 9 2
H           -3.884864        2.369505        3.794311 9 3
H            0.804352       -6.170201       -2.455154 10 2
H            0.489826       -5.114382       -1.417222 10 3
H            1.096538       -1.240000      -10.685961 11 2
H            0.070907       -1.650526       -9.651296 11 3
H           -7.062196       -1.010451       -2.708807 12 2
H           -8.106394       -1.340645       -1.664013 12 3
%EOF%
#######################################################################
#  QM/MM INTERACTION INPUT
#######################################################################
cat > ccsdmm_modelmix.pot <<%EOF%
**SYSTP
.NUMMTP   # It is possible to mix the different models. As SPC models
 6        # does not couple the t and t-bar parameters they are much
.TYPE     # faster than the full SPC_E01 model. As the polarization
 0        # interaction is proportional to 1/R^3 they fall off much
.MODEL    # faster than the coulomb interaction which is proportional
 SPC_E01  # to 1/R. It might be enough to include a polarization as
.CHARGS   # SPC_E01 only on the closest solvent molecules and SPC(_EC#)
 3        # on the rest. The scaling of the SPC_E01 iterations are
 -0.669   # N^5 (like CC2) where N is the number of basis functions.
 0.3345
 0.3345
.ALPISO
 3
 9.500   # It is possible to give atomic polarizabilities. Here, 
 0.200   # alpha_O=9.500, alpha_H=0.200
 0.200   # The order is the same as the sub system site order. See MOLECULE.INP!
*******
.TYPE
 1-2
.MODEL
 SPC_E01
.ALPISO
 1
 9.718
*******
.TYPE
 3-4
.MODEL
 SPC_E01
.ALPISO
 3
 5.340
 2.300
 2.300
*******
.TYPE
 5-6
.MODEL
 SPC_EC3
.ALPISO
 1
 9.718
*******
.TYPE
 7-8
.MODEL
 SPC_EC3
.ALPISO
 3
 5.340
 2.300
 2.300
*******
.TYPE
 9-10
.MODEL
 SPC_EC1
.ALPISO
 3
 5.340
 2.300
 2.300
*******
.TYPE
 11-12
.MODEL
 SPC
*******
**TWOIA (i,j=0,1,2,...,N; if i=0 then j.neq.0)
.LJ_A
 27           # System 1     System 2
 2083000      #        0            1
 2083000      #        0            2
 2083000      #        0            3
 2083000      #        0            4
 2083000      #        0            5
 2083000      #        0            6
 2083000      #        1            1
 2083000      #        1            2
 2083000      #        1            3
 2083000      #        1            4
 2083000      #        1            5
 2083000      #        1            6
 2083000      #        2            2
 2083000      #        2            3
 2083000      #        2            4
 2083000      #        2            5
 2083000      #        2            6
 2083000      #        3            3
 2083000      #        3            4
 2083000      #        3            5
 2083000      #        3            6
 2083000      #        4            4
 2083000      #        4            5
 2083000      #        4            6
 2083000      #        5            5
 2083000      #        5            6
 2083000      #        6            6  And the same for B parameters
.LJ_B
 27           # The total number of parameters:
 45.21        # 2*(# MMTYPES) + (MMTYPES - 1) + (MMTYPES - 2) + 
 45.21        # (MMTYPES - 3) + ... + 1
 45.21        # MMTYPES=6 as in this case gives
 45.21        # 2*6 + 5 + 4 + 3 + 2 + 1 = 27
 45.21
 45.21
 45.21
 45.21
 45.21
 45.21
 45.21
 45.21
 45.21
 45.21
 45.21
 45.21
 45.21
 45.21
 45.21
 45.21
 45.21
 45.21
 45.21 
 45.21
 45.21
 45.21
 45.21
**END OF
%EOF%
#
#######################################################################
#  DALTON INPUT
#######################################################################
cat > ccsdmm_modelmix.dal <<%EOF%
**DALTON INPUT
.RUN WAVEFUNCTION
*QM3
.QM3
.THRDIP
 1.0D-11
.MAXDIP
 80
.OLDTG
**INTEGRALS
.DIPLEN
.NUCPOT
.NELFLD
.THETA
.SECMOM
**WAVE FUNCTIONS
.CC
.PRINT
 6
*SCF INPUT
.THRESH
1.0D-11
*CC INP
.CCSD
.THRLEQ
 1.0D-11
.THRENR
 1.0D-11
.MAX IT
 90
.MXLRV
 180
*CCSLV
.CCMM
.ETOLSL
 1.0D-10
.TTOLSL
 1.0D-10
.LTOLSL
 1.0D-10
.MXSLIT
 200
.MXINIT
 4 5
*CCFOP
.DIPMOM
.QUADRU
.SECMOM
.NONREL
*CCLR
.DIPOLE
.OLD_LR
.FREQUE
 1
 0.0000
*CCEXCI
.NCCEXCI
 1 0 0 0
*CCLRSD
.DIPOLE
**END OF
%EOF%
#######################################################################
#  CHECK SCRIPT
#######################################################################
echo $CHECK_SHELL >ccsdmm_modelmix.check
cat >>ccsdmm_modelmix.check <<'%EOF%'
log=$1

if [ `uname` = Linux ]; then
   GREP="egrep -a"
else
   GREP="egrep"
fi

if $GREP -q "not implemented for parallel calculations" $log; then
   echo "TEST ENDED AS EXPECTED"
   exit 0
fi

# QM/MM interaction energy compared:
CRIT1=`$GREP "\| * (\-|\-0).010193083. \| * (\-|\-0)\.001413232. \| * ( |0)\.005071554. \| * (\-|\-0)\.00653493.. \|" $log | wc -l`
CRIT2=`$GREP "\| * \-75\.71407834.. \| * \-75\.72061327.. \| * (\-|\-0)\.000001850. \|" $log | wc -l`
TEST[1]=`expr $CRIT1 \+ $CRIT2`
CTRL[1]=2
ERROR[1]="QM/MM ENERGY NOT CORRECT"

# Dipole moment components compared:
CRIT1=`$GREP "x * ( |0)\.09378868 * ( |0)\.23838704" $log | wc -l`
CRIT2=`$GREP "y * ( |0)\.00318384 * ( |0)\.00809251" $log | wc -l`
CRIT3=`$GREP "z * ( |0)\.94747237 * 2\.40823439" $log | wc -l` 
TEST[2]=`expr $CRIT1 \+ $CRIT2 \+ $CRIT3`
CTRL[2]=7
ERROR[2]="DIPOLE MOMENT NOT CORRECT"

# Quadrupole moment components compared:
CRIT1=`$GREP "1 * 1\.3936093. * ( |0)\.0004815. * ( |0)\.1478202." $log | wc -l`
CRIT2=`$GREP "2 *  ( |0)\.0004815.  *  -1\.5681122. * ( |0)\.0014056." $log | wc -l`
CRIT3=`$GREP "3 * ( |0)\.1478202. * ( |0)\.0014056. * ( |0)\.1745028." $log | wc -l`
TEST[3]=`expr $CRIT1 \+ $CRIT2 \+ $CRIT3`
CTRL[3]=9
ERROR[3]="QUADRUPOLE MOMENT NOT CORRECT"

# Second order electric moment components compared:
CRIT1=`$GREP "1 * 7\.1879348. * (\-|\-0)\.00032101 * (\-|\-0)\.09854685" $log | wc -l`
CRIT2=`$GREP "2 * (\-|\-0)\.00032101 * 5\.0717988. * (\-|\-0)\.00093713" $log | wc -l`
CRIT3=`$GREP "3 * (\-|\-0)\.09854685 * (\-|\-0)\.00093713 * 6\.3626927." $log | wc -l`
TEST[4]=`expr $CRIT1 \+ $CRIT2 \+ $CRIT3`
CTRL[4]=10
ERROR[4]="SECOND ORDER MOMENT NOT CORRECT"

# First ground state excitation energy compared:
CRIT1=`$GREP "1A * \| * 1 * \| * ( |0)\.3265... * \| *  8\.8849. * \| *  71662\...." $log | wc -l`
TEST[5]=`expr $CRIT1`
CTRL[5]=1
ERROR[5]="FIRST GROUND STATE EXCITATION ENERGY NOT CORRECT"

# Transition moments between ground state and first excited state compared:
CRIT1=`$GREP "1A * \| * 1 * \| * ( |0)\.0265... * \| * ( |0)\.0057... * \|" $log | wc -l`
TEST[6]=`expr $CRIT1`
CTRL[6]=1
ERROR[6]="OSCILLATOR STRENGTH NOT CORRECT"

# Static polarization compared:
CRIT1=`$GREP "1 *  6\.355531.. * (\-|\-0)\.000866.. * (\-|\-0)\.256175.." $log | wc -l`
CRIT2=`$GREP "2 * (\-|\-0)\.000866.. * ( |0)\.818662.. * (\-|\-0)\.004419.." $log | wc -l` 
CRIT3=`$GREP "3 * (\-|\-0)\.256175.. * (\-|\-0)\.004419.. * 3\.781498.." $log | wc -l`
TEST[7]=`expr $CRIT1 \+ $CRIT2 \+ $CRIT3`
CTRL[7]=6
ERROR[7]="STATIC POLARIZATION NOT CORRECT"

PASSED=1
for i in 1 2 3 4 5 6 7
do 
   if [ ${TEST[i]} -ne ${CTRL[i]} ]; then
     echo "${ERROR[i]} ( test = ${TEST[i]}; control = ${CTRL[i]} ); "
     PASSED=0
   fi
done 

if [ $PASSED -eq 1 ]
then
  echo TEST ENDED PROPERLY
  exit 0
else
  echo THERE IS A PROBLEM 
  exit 1
fi
%EOF%
#######################################################################
